library(data.table)
library(ggplot2)
library(here)
library(gridExtra)
library(grid)
library(gtable)
dir = here()
wd = gsub("figures", "data/world bank", dir)
setwd(wd)
d = fread("data_format.csv")
d$EG.USE.PCAP.KG.OE = d$EG.USE.PCAP.KG.OE*4.187/100
text.size = 10
# intensity
#############################################################
# regression
x = log(d$SL.SRV.EMPL.ZS)
y = log(d$EG.EGY.PRIM.PP.KD)
x = x[d$EG.USE.COMM.FO.ZS > 0]
y = y[d$EG.USE.COMM.FO.ZS > 0]
r = lm(y ~ x)
r.square = summary(r)$r.squared
x = seq(log(1), log(100), length.out = 50)
x = data.frame(x)
regression = exp(predict(r, x,  interval = "prediction", level = 0.99))
regression = data.frame(x = exp(x), regression)
intensity.r2 =  paste("R^2 == ", round(r.square,2))
intensity = ggplot() +
geom_ribbon(data = regression, aes( x = x, ymin = lwr, ymax = upr), fill = "black", alpha = 0.05) +
geom_line(data = regression, aes(x = x, y = fit), col = "grey50") +
geom_point(data = d, aes(x = SL.SRV.EMPL.ZS, y = EG.EGY.PRIM.PP.KD), size = 0.1, col = "dodgerblue4") +
scale_x_log10("Service Sector (% of Employment)", breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100)) +
scale_y_log10("MJ per $2011 PPP GDP", breaks = c(0.1, 1, 10, 100), labels = c(0.1, 1, 10, 100)) +
coord_cartesian(xlim = c(4, 100), ylim = c(0.1, 100)) +
ggtitle("A. Energy Intensity of GDP") +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5, vjust = -0.1),
axis.line = element_line(color = "black"),
axis.title.x=element_text(size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times")) +
annotate("text", x = 5, y = 100, label = intensity.r2, parse = T,family="Times", size=3)
# fossil fuel fraction
#############################################################
# regression
x = log(d$SL.SRV.EMPL.ZS)
y = log(d$EG.USE.COMM.FO.ZS)
x = x[d$EG.USE.COMM.FO.ZS > 0]
y = y[d$EG.USE.COMM.FO.ZS > 0]
r = lm(y ~ x)
r.square = summary(r)$r.squared
x = seq(log(1), log(100), length.out = 50)
x = data.frame(x)
regression = exp(predict(r, x,  interval = "prediction", level = 0.99))
regression = data.frame(x = exp(x), regression)
fossil.r2 =  paste("R^2 == ", round(r.square,2))
fossil = ggplot() +
geom_ribbon(data = regression, aes( x = x, ymin = lwr, ymax = upr), fill = "black", alpha = 0.05) +
geom_line(data = regression, aes(x = x, y = fit), col = "grey50") +
geom_point(data = d, aes(x = SL.SRV.EMPL.ZS, y = EG.USE.COMM.FO.ZS), size = 0.1, col = "dodgerblue4") +
scale_x_log10("Service Sector (% of Employment)", breaks = c(2, 5, 10, 20, 50, 100)) +
scale_y_log10("Fossil Fuels (% Total Energy)",
breaks = c(2, 5, 10, 20, 50, 100) ) +
coord_cartesian(xlim = c(4, 100), ylim = c(1.5, 100)) +
ggtitle("B.  Percentage of Energy from Fossil Fuels") +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5, vjust = -0.1),
axis.line = element_line(color = "black"),
axis.title.x=element_text(size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times")) +
annotate("text", x = 5, y = 100, label = fossil.r2, parse = T,family="Times", size=3)
# total energy
####################################################################
# regression
x = log(d$SL.SRV.EMPL.ZS)
y = log(d$EG.USE.PCAP.KG.OE)
r = lm(y ~ x)
r.square = summary(r)$r.squared
x = seq(log(1), log(100), length.out = 50)
x = data.frame(x)
regression = exp(predict(r, x,  interval = "prediction", level = 0.99))
regression = data.frame(x = exp(x), regression)
energy.r2 =  paste("R^2 == ", round(r.square,2))
energy = ggplot() +
geom_ribbon(data = regression, aes( x = x, ymin = lwr, ymax = upr), fill = "black", alpha = 0.05) +
geom_line(data = regression, aes(x = x, y = fit), col = "grey50") +
geom_point(data = d, aes(x = SL.SRV.EMPL.ZS, y = EG.USE.PCAP.KG.OE), size = 0.1, col = "dodgerblue4") +
scale_x_log10("Service Sector (% of Employment)", breaks = c(2, 5, 10, 20, 50, 100)) +
scale_y_log10("Energy Use per Capita (GJ)", breaks =  c(1,10, 100, 1000)  ) +
coord_cartesian(xlim = c(4, 100), ylim = c(0.4, 1000)) +
ggtitle("C.  Energy Use per Capita") +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5, vjust = -0.1),
axis.line = element_line(color = "black"),
axis.title.x=element_text(size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times")) +
annotate("text", x = 5, y = 1000, label = energy.r2, parse = T,family="Times", size=3)
# global energy
######################################################################
g =  fread("sector_energy.csv")
wd
